//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : CudaWavesComputer.cu
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#include <math.h>
#include <cuda.h>
#include <cublas.h>
#include <stdio.h>
#include <sys/time.h>
#include "assertplus.h"
#include "FloatUtils.h"
#include "CudaWavesComputer.h"

#include "CUDAPACKAGE/common/defs.h"
#include "CUDAPACKAGE/common/macros.h"
#include "CUDAPACKAGE/common/misc.h"

// Time integration / waves specific CUDA code
#include "CUDAPACKAGE/wavescomputer/reorder.h"
#include "CUDAPACKAGE/wavescomputer/shiftzetaphi.h"
#include "CUDAPACKAGE/wavescomputer/boundaryzetaphi.h"
#include "CUDAPACKAGE/wavescomputer/solvepsi.h"
#include "CUDAPACKAGE/wavescomputer/adjustpsi.h"
#include "CUDAPACKAGE/wavescomputer/shiftpsi.h"
#include "CUDAPACKAGE/wavescomputer/wavebreaking.h"
#include "CUDAPACKAGE/wavescomputer/boundarypsi.h"
#include "CUDAPACKAGE/wavescomputer/timederivatives.h"
#include "CUDAPACKAGE/wavescomputer/smoothzetaphi.h"
#include "CUDAPACKAGE/wavescomputer/leapfrog.h"


/// constructor
CudaWavesComputer::CudaWavesComputer(const int Nx1,
                                     const int Nx2,
                                     const REAL dx1,
                                     const REAL dx2,
                                     const int bcN,
                                     const int bcE,
                                     const int bcS,
                                     const int bcW,
                                     const REAL g)
: m_Nx1(Nx1)
, m_Nx2(Nx2)
, m_dx1(dx1)
, m_dx2(dx2)
, m_bcN(bcN)
, m_bcE(bcE)
, m_bcS(bcS)
, m_bcW(bcW)
, m_g(g)
, m_Nc(0)
, m_Nr(0)
, m_n(0)
, m_dpsi(0)
, m_dpsim1(0)
, m_dpsim2(0)
, m_dpsi_in_N_ext(0)
, m_dpsi_in_N_cen(0)
, m_dpsiN(0)
, m_dpsi_in_E_ext(0)
, m_dpsi_in_E_cen(0)
, m_dpsiE(0)
, m_dpsi_in_S_ext(0)
, m_dpsi_in_S_cen(0)
, m_dpsiS(0)
, m_dpsi_in_W_ext(0)
, m_dpsi_in_W_cen(0)
, m_dpsiW(0)
, m_dphi_in_N_ext(0)
, m_dphi_in_N_cen(0)
, m_dphiN(0)
, m_dphi_in_E_ext(0)
, m_dphi_in_E_cen(0)
, m_dphiE(0)
, m_dphi_in_S_ext(0)
, m_dphi_in_S_cen(0)
, m_dphiS(0)
, m_dphi_in_W_ext(0)
, m_dphi_in_W_cen(0)
, m_dphiW(0)
, m_dzeta_in_N_ext(0)
, m_dzeta_in_N_cen(0)
, m_dzetaN(0)
, m_dzeta_in_E_ext(0)
, m_dzeta_in_E_cen(0)
, m_dzetaE(0)
, m_dzeta_in_S_ext(0)
, m_dzeta_in_S_cen(0)
, m_dzetaS(0)
, m_dzeta_in_W_ext(0)
, m_dzeta_in_W_cen(0)
, m_dzetaW(0)
, m_dstoN(0)
, m_dstoE(0)
, m_dstoS(0)
, m_dstoW(0)
, m_dbufc(0) // for reordering of .. x 3 arrays
, m_dbufr(0)
, m_dpressurePulse(0)
, m_dU1(0)
, m_dU2(0)
, m_dh0(0)
, m_dwet(0)
, m_dD0(0)
, m_dzeta(0)
, m_dzetam1(0)
, m_dzetam2(0)
, m_ddt_zeta(0) // derivative d zeta / dt
, m_ddt_zetam1(0)
, m_dphi(0)
, m_dphim1(0)
, m_dphim2(0)
, m_ddt_phi(0)  // derivative d phi / dt
, m_ddt_phim1(0)
, m_dSRCz(0)
, m_dSRCp(0) // used in solvePsi()
, m_dbpsi(0) // RHS in linear system S*psi = b
, m_dSppC(0) // <-- what is this?
, m_dalfa_brk(0)
, m_dzeta_std(0)
, m_numBlocks(0)
, m_numThreads(256) // number of CUDA threads
{
    determineSizes();

    setup();
} // constructor


/// destructor
CudaWavesComputer::~CudaWavesComputer()
{
    cleanup();
} // destructor


void CudaWavesComputer::determineSizes()
{
    m_nx[0] = DIM_COMPUTE_BLOCK * (int)ceil((float)m_Nx2 / DIM_COMPUTE_BLOCK) 
               + BORDER_WIDTH2;

    m_ny[0] = DIM_COMPUTE_BLOCK * (int)ceil((float)m_Nx1 / DIM_COMPUTE_BLOCK) 
               + BORDER_WIDTH2;
  
    m_n = m_nx[0] * m_ny[0];

    m_cx[0] = m_nx[0] - BORDER_WIDTH2;
    m_cy[0] = m_ny[0] - BORDER_WIDTH2;

    m_Nc = NUM_THREADS * (int)ceil((float)(m_Nx1) / NUM_THREADS); // padded
    m_Nr = NUM_THREADS * (int)ceil((float)(m_Nx2) / NUM_THREADS); // padded
} // determineSizes


void CudaWavesComputer::setup()
{
    // Array2D(COLUMNS, ROWS) objects (same embedding format as CudaRrbSolver)
    m_dpsi           = common::initGpuVector<REAL,1>(m_n, 0.0); // if flag = 1: initialize with 0.0
    m_dpsim1         = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dpsim2         = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dU1            = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dU2            = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dh0            = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dwet           = common::initGpuVector<SHORT,1>(m_n, 0.0);
    m_dD0            = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dzeta          = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_ddt_zeta       = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_ddt_zetam1     = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dphi           = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_ddt_phi        = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_ddt_phim1      = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dSRCz          = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dpressurePulse = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dzetam1        = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dzetam2        = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dphim1         = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dphim2         = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dalfa_brk      = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dzeta_std      = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dSRCp          = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dbpsi          = common::initGpuVector<REAL,1>(m_n, 0.0);
    m_dSppC          = common::initGpuVector<REAL,1>(m_n, 0.0);

    // Array2D(COLUMNS, 3) objects & Array2D(ROWS, 3) objects
    m_dbufc          = common::initGpuVector<REAL,1>(m_Nc * 3, 0.0); // buffer for reordering column vectors
    m_dbufr          = common::initGpuVector<REAL,1>(m_Nr * 3, 0.0); // buffer for reordering row vector
    m_dpsiN          = common::initGpuVector<REAL,1>(m_Nc * 3, 0.0);
    m_dpsiE          = common::initGpuVector<REAL,1>(m_Nr * 3, 0.0);
    m_dpsiS          = common::initGpuVector<REAL,1>(m_Nc * 3, 0.0);
    m_dpsiW          = common::initGpuVector<REAL,1>(m_Nr * 3, 0.0);
    m_dphiN          = common::initGpuVector<REAL,1>(m_Nc * 3, 0.0);
    m_dphiE          = common::initGpuVector<REAL,1>(m_Nr * 3, 0.0);
    m_dphiS          = common::initGpuVector<REAL,1>(m_Nc * 3, 0.0);
    m_dphiW          = common::initGpuVector<REAL,1>(m_Nr * 3, 0.0);
    m_dzetaN         = common::initGpuVector<REAL,1>(m_Nc * 3, 0.0);
    m_dzetaE         = common::initGpuVector<REAL,1>(m_Nr * 3, 0.0);
    m_dzetaS         = common::initGpuVector<REAL,1>(m_Nc * 3, 0.0);
    m_dzetaW         = common::initGpuVector<REAL,1>(m_Nr * 3, 0.0);

    // Array1D(COLUMNS) & Array1D(ROWS) objects
    m_dpsi_in_N_ext  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dpsi_in_N_cen  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dpsi_in_E_ext  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dpsi_in_E_cen  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dpsi_in_S_ext  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dpsi_in_S_cen  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dpsi_in_W_ext  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dpsi_in_W_cen  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dphi_in_N_ext  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dphi_in_N_cen  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dphi_in_E_ext  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dphi_in_E_cen  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dphi_in_S_ext  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dphi_in_S_cen  = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dphi_in_W_ext  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dphi_in_W_cen  = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dzeta_in_N_ext = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dzeta_in_N_cen = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dzeta_in_E_ext = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dzeta_in_E_cen = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dzeta_in_S_ext = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dzeta_in_S_cen = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dzeta_in_W_ext = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dzeta_in_W_cen = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dstoN          = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dstoE          = common::initGpuVector<REAL,1>(m_Nr, 0.0);
    m_dstoS          = common::initGpuVector<REAL,1>(m_Nc, 0.0);
    m_dstoW          = common::initGpuVector<REAL,1>(m_Nr, 0.0);
} // setup


void CudaWavesComputer::cleanup()
{
    // Array2D(COLUMNS, ROWS) objects (same embedding format as CudaRrbSolver)
    common::freeGpuVector<REAL>(m_dpsi);
    common::freeGpuVector<REAL>(m_dpsim1);
    common::freeGpuVector<REAL>(m_dpsim2);
    common::freeGpuVector<REAL>(m_dU1);
    common::freeGpuVector<REAL>(m_dU2);
    common::freeGpuVector<REAL>(m_dh0);
    common::freeGpuVector<SHORT>(m_dwet);
    common::freeGpuVector<REAL>(m_dD0);
    common::freeGpuVector<REAL>(m_dzeta);
    common::freeGpuVector<REAL>(m_ddt_zeta);
    common::freeGpuVector<REAL>(m_ddt_zetam1);
    common::freeGpuVector<REAL>(m_dphi);
    common::freeGpuVector<REAL>(m_ddt_phi);
    common::freeGpuVector<REAL>(m_ddt_phim1);
    common::freeGpuVector<REAL>(m_dSRCz);
    common::freeGpuVector<REAL>(m_dpressurePulse);
    common::freeGpuVector<REAL>(m_dzetam1);
    common::freeGpuVector<REAL>(m_dzetam2);
    common::freeGpuVector<REAL>(m_dphim1);
    common::freeGpuVector<REAL>(m_dphim2);
    common::freeGpuVector<REAL>(m_dalfa_brk);
    common::freeGpuVector<REAL>(m_dzeta_std);
    common::freeGpuVector<REAL>(m_dSRCp);
    common::freeGpuVector<REAL>(m_dbpsi);
    common::freeGpuVector<REAL>(m_dSppC);

    // Array2D(COLUMNS, 3) objects & Array2D(ROWS, 3) objects
    common::freeGpuVector<REAL>(m_dbufc);
    common::freeGpuVector<REAL>(m_dbufr);
    common::freeGpuVector<REAL>(m_dpsiN);
    common::freeGpuVector<REAL>(m_dpsiE);
    common::freeGpuVector<REAL>(m_dpsiS);
    common::freeGpuVector<REAL>(m_dpsiW);
    common::freeGpuVector<REAL>(m_dphiN);
    common::freeGpuVector<REAL>(m_dphiE);
    common::freeGpuVector<REAL>(m_dphiS);
    common::freeGpuVector<REAL>(m_dphiW);
    common::freeGpuVector<REAL>(m_dzetaN);
    common::freeGpuVector<REAL>(m_dzetaE);
    common::freeGpuVector<REAL>(m_dzetaS);
    common::freeGpuVector<REAL>(m_dzetaW);

    // Array1D(COLUMNS) & Array1D(ROWS) objects
    common::freeGpuVector<REAL>(m_dpsi_in_N_ext);
    common::freeGpuVector<REAL>(m_dpsi_in_N_cen);
    common::freeGpuVector<REAL>(m_dpsi_in_E_ext);
    common::freeGpuVector<REAL>(m_dpsi_in_E_cen);
    common::freeGpuVector<REAL>(m_dpsi_in_S_ext);
    common::freeGpuVector<REAL>(m_dpsi_in_S_cen);
    common::freeGpuVector<REAL>(m_dpsi_in_W_ext);
    common::freeGpuVector<REAL>(m_dpsi_in_W_cen);
    common::freeGpuVector<REAL>(m_dphi_in_N_ext);
    common::freeGpuVector<REAL>(m_dphi_in_N_cen);
    common::freeGpuVector<REAL>(m_dphi_in_E_ext);
    common::freeGpuVector<REAL>(m_dphi_in_E_cen);
    common::freeGpuVector<REAL>(m_dphi_in_S_ext);
    common::freeGpuVector<REAL>(m_dphi_in_S_cen);
    common::freeGpuVector<REAL>(m_dphi_in_W_ext);
    common::freeGpuVector<REAL>(m_dphi_in_W_cen);
    common::freeGpuVector<REAL>(m_dzeta_in_N_ext);
    common::freeGpuVector<REAL>(m_dzeta_in_N_cen);
    common::freeGpuVector<REAL>(m_dzeta_in_E_ext);
    common::freeGpuVector<REAL>(m_dzeta_in_E_cen);
    common::freeGpuVector<REAL>(m_dzeta_in_S_ext);
    common::freeGpuVector<REAL>(m_dzeta_in_S_cen);
    common::freeGpuVector<REAL>(m_dzeta_in_W_ext);
    common::freeGpuVector<REAL>(m_dzeta_in_W_cen);
    common::freeGpuVector<REAL>(m_dstoN);
    common::freeGpuVector<REAL>(m_dstoE);
    common::freeGpuVector<REAL>(m_dstoS);
    common::freeGpuVector<REAL>(m_dstoW);
} // cleanup


// arrays that are constant during simulation are copied during initialization
// of the wavescomputer to the device
void CudaWavesComputer::initialize(const Array2D<REAL> *hU1,
                                   const Array2D<REAL> *hU2,
                                   const Array2D<REAL> *hh0,
                                   const Array2D<SHORT> *hwet,
                                   const Array2D<REAL> *hD0,
                                   const Array2D<REAL> *hSRCz,
                                   const Array2D<REAL> *hSRCp,
                                   const Array2D<REAL> *halfa_brk,
                                   const Array2D<REAL> *hzeta_std,
                                   const Array2D<REAL> *hSppC,
                                   const Array1D<REAL> *hstoN,
                                   const Array1D<REAL> *hstoE,
                                   const Array1D<REAL> *hstoS,
                                   const Array1D<REAL> *hstoW)
{
    // Array2D objects:
    common::setGpu2DVector<REAL>(m_dU1, hU1, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dU2, hU2, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dh0, hh0, m_nx[0]);
    common::setGpu2DVector<SHORT>(m_dwet, hwet, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dD0, hD0, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dSRCz, hSRCz, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dSRCp, hSRCp, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dalfa_brk, halfa_brk, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dzeta_std, hzeta_std, m_nx[0]);
    common::setGpu2DVector<REAL>(m_dSppC, hSppC, m_nx[0]);

    // Array1D objects:
    common::setGpu1DVector<REAL>(m_dstoN, hstoN);
    common::setGpu1DVector<REAL>(m_dstoE, hstoE);
    common::setGpu1DVector<REAL>(m_dstoS, hstoS);
    common::setGpu1DVector<REAL>(m_dstoW, hstoW);
} // intialize


void CudaWavesComputer::shiftZetaPhi(const REAL eps)
{
        wavescomputer::shiftzetaphi<REAL>
           (m_dzetam2, m_dzetam1, m_dzeta,
            m_dphim2, m_dphim1, m_dphi,
            m_ddt_zetam1, m_ddt_zeta,
            m_ddt_phim1, m_ddt_phi,
            eps,
            m_nx[0], m_ny[0]);
}


void CudaWavesComputer::waveBreaking(Array2D<REAL> *hzetastd,
                                     Array2D<REAL> *halfabrk,
                                     REAL zeta_mean, REAL tau_brk, REAL dt_new)
{
    wavescomputer::wavebreaking<REAL>(m_dzeta_std,
                                      m_dalfa_brk,
                                      m_dzeta,
                                      m_dwet,
                                      m_dh0,
                                      zeta_mean,
                                      tau_brk,
                                      dt_new,
                                      m_nx[0],
                                      m_ny[0]);
    // nagaan of dit wel nodig is:
    common::getGpu2DVector<REAL>(halfabrk,  m_dalfa_brk,  m_nx[0]);
    common::getGpu2DVector<REAL>(hzetastd,  m_dzeta_std,  m_nx[0]);
}

void CudaWavesComputer::boundaryZetaPhi(Array2D<REAL> *hphi,
                                        Array2D<REAL> *hzeta,
                                        const Array1D<REAL> *hphi_in_N_ext,
                                        const Array1D<REAL> *hphi_in_N_cen,
                                        const Array2D<REAL> *hphiN,
                                        const Array1D<REAL> *hzeta_in_N_ext,
                                        const Array1D<REAL> *hzeta_in_N_cen,
                                        const Array2D<REAL> *hzetaN,
                                        const Array1D<REAL> *hphi_in_E_ext,
                                        const Array1D<REAL> *hphi_in_E_cen,
                                        const Array2D<REAL> *hphiE,
                                        const Array1D<REAL> *hzeta_in_E_ext,
                                        const Array1D<REAL> *hzeta_in_E_cen,
                                        const Array2D<REAL> *hzetaE,
                                        const Array1D<REAL> *hphi_in_S_ext,
                                        const Array1D<REAL> *hphi_in_S_cen,
                                        const Array2D<REAL> *hphiS,
                                        const Array1D<REAL> *hzeta_in_S_ext,
                                        const Array1D<REAL> *hzeta_in_S_cen,
                                        const Array2D<REAL> *hzetaS,
                                        const Array1D<REAL> *hphi_in_W_ext,
                                        const Array1D<REAL> *hphi_in_W_cen,
                                        const Array2D<REAL> *hphiW,
                                        const Array1D<REAL> *hzeta_in_W_ext,
                                        const Array1D<REAL> *hzeta_in_W_cen,
                                        const Array2D<REAL> *hzetaW)
{
    // North
    common::setGpuBdVector<REAL>(m_dbufc, hphiN);
    wavescomputer::reorder<REAL>(m_dphiN,  m_dbufc, m_Nx1);

    common::setGpuBdVector<REAL>(m_dbufc, hzetaN);
    wavescomputer::reorder<REAL>(m_dzetaN, m_dbufc, m_Nx1);


    // East
    common::setGpuBdVector<REAL>(m_dbufr, hphiE);
    wavescomputer::reorder<REAL>(m_dphiE,  m_dbufr, m_Nx2);

    common::setGpuBdVector<REAL>(m_dbufr, hzetaE);
    wavescomputer::reorder<REAL>(m_dzetaE, m_dbufr, m_Nx2);


    // South
    common::setGpuBdVector<REAL>(m_dbufc, hphiS);
    wavescomputer::reorder<REAL>(m_dphiS,  m_dbufc, m_Nx1);

    common::setGpuBdVector<REAL>(m_dbufc, hzetaS);
    wavescomputer::reorder<REAL>(m_dzetaS, m_dbufc, m_Nx1);


    // West
    common::setGpuBdVector<REAL>(m_dbufr, hphiW);;
    wavescomputer::reorder<REAL>(m_dphiW,  m_dbufr, m_Nx2);

    common::setGpuBdVector<REAL>(m_dbufr, hzetaW);
    wavescomputer::reorder<REAL>(m_dzetaW, m_dbufr, m_Nx2);


    // North
    if (m_bcN == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_N_ext,  hphi_in_N_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_N_ext, hzeta_in_N_ext);
        common::setGpu1DVector<REAL>(m_dphi_in_N_cen,  hphi_in_N_cen);
        common::setGpu1DVector<REAL>(m_dzeta_in_N_cen, hzeta_in_N_cen);

        wavescomputer::boundaryzetaphi::reflectiveneumann::north<REAL>
            (m_dphi, m_dzeta,
             m_dphiN, m_dzetaN, 
             m_dphi_in_N_ext, m_dzeta_in_N_ext,
             m_dphi_in_N_cen, m_dzeta_in_N_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcN == 1) // weakly-reflective boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_N_ext,  hphi_in_N_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_N_ext, hzeta_in_N_ext);

        wavescomputer::boundaryzetaphi::weaklyreflective::north<REAL>
            (m_dphi, m_dzeta,
             m_dphiN, m_dzetaN,
             m_dphi_in_N_ext, m_dzeta_in_N_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundaryzetaphi<north>: unknown boundary condition");
    }


    // East
    if (m_bcE == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_E_ext,  hphi_in_E_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_E_ext, hzeta_in_E_ext);
        common::setGpu1DVector<REAL>(m_dphi_in_E_cen,  hphi_in_E_cen);
        common::setGpu1DVector<REAL>(m_dzeta_in_E_cen, hzeta_in_E_cen);

        wavescomputer::boundaryzetaphi::reflectiveneumann::east<REAL>
            (m_dphi, m_dzeta,
             m_dphiE, m_dzetaE, 
             m_dphi_in_E_ext, m_dzeta_in_E_ext,
             m_dphi_in_E_cen, m_dzeta_in_E_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcE == 1) // weakly-reflective boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_E_ext,  hphi_in_E_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_E_ext, hzeta_in_E_ext);

        wavescomputer::boundaryzetaphi::weaklyreflective::east<REAL>
            (m_dphi, m_dzeta,
             m_dphiE, m_dzetaE,
             m_dphi_in_E_ext, m_dzeta_in_E_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundaryzetaphi<east>: unknown boundary condition");
    }


    // South
    if (m_bcS == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_S_ext, hphi_in_S_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_S_ext, hzeta_in_S_ext);
        common::setGpu1DVector<REAL>(m_dphi_in_S_cen, hphi_in_S_cen);
        common::setGpu1DVector<REAL>(m_dzeta_in_S_cen, hzeta_in_S_cen);

        wavescomputer::boundaryzetaphi::reflectiveneumann::south<REAL>
            (m_dphi, m_dzeta,
             m_dphiS, m_dzetaS, 
             m_dphi_in_S_ext, m_dzeta_in_S_ext,
             m_dphi_in_S_cen, m_dzeta_in_S_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcS == 1) // weakly-reflective boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_S_ext, hphi_in_S_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_S_ext, hzeta_in_S_ext);
        
        wavescomputer::boundaryzetaphi::weaklyreflective::south<REAL>
            (m_dphi, m_dzeta,
             m_dphiS, m_dzetaS,
             m_dphi_in_S_ext, m_dzeta_in_S_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundaryzetaphi<south>: unknown boundary condition");
    }


    // West
    if (m_bcW == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_W_ext,  hphi_in_W_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_W_ext, hzeta_in_W_ext);
        common::setGpu1DVector<REAL>(m_dphi_in_W_cen,  hphi_in_W_cen);
        common::setGpu1DVector<REAL>(m_dzeta_in_W_cen, hzeta_in_W_cen);
  
        wavescomputer::boundaryzetaphi::reflectiveneumann::west<REAL>
            (m_dphi, m_dzeta,
             m_dphiW, m_dzetaW, 
             m_dphi_in_W_ext, m_dzeta_in_W_ext,
             m_dphi_in_W_cen, m_dzeta_in_W_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcW == 1) // weakly-reflective boundary
    {
        common::setGpu1DVector<REAL>(m_dphi_in_W_ext,  hphi_in_W_ext);
        common::setGpu1DVector<REAL>(m_dzeta_in_W_ext, hzeta_in_W_ext);
         
        wavescomputer::boundaryzetaphi::weaklyreflective::west<REAL>
            (m_dphi, m_dzeta,
             m_dphiW, m_dzetaW,
             m_dphi_in_W_ext, m_dzeta_in_W_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundaryzetaphi<west>: unknown boundary condition");
    }

    common::getGpu2DVector<REAL>(hphi,  m_dphi,  m_nx[0]); // needed for weaklyReflective()
    common::getGpu2DVector<REAL>(hzeta, m_dzeta, m_nx[0]); // needed for weaklyReflective()
}


void CudaWavesComputer::solvePsi(Array2D<REAL> *hbpsi,
                                 const REAL hpslow,
                                 const Array2D<REAL> *hpsiN,
                                 const Array1D<REAL> *hpsi_in_N_ext,
                                 const Array1D<REAL> *hpsi_in_N_cen,
                                 const Array2D<REAL> *hpsiE,
                                 const Array1D<REAL> *hpsi_in_E_ext,
                                 const Array1D<REAL> *hpsi_in_E_cen,
                                 const Array2D<REAL> *hpsiS,
                                 const Array1D<REAL> *hpsi_in_S_ext,
                                 const Array1D<REAL> *hpsi_in_S_cen,
                                 const Array2D<REAL> *hpsiW,
                                 const Array1D<REAL> *hpsi_in_W_ext,
                                 const Array1D<REAL> *hpsi_in_W_cen)
{
    // Step 1) compute bpsi for entire field
    wavescomputer::solvepsi::phase1<REAL>(m_dbpsi, m_dh0, m_dD0, 
                                  m_dphi, m_dSRCp, m_dSppC, 
                                  m_nx[0], m_Nx1, m_Nx2,
                                  m_cx[0], m_cy[0],
                                  m_dx1, m_dx2, 
                                  hpslow);


    // Step 2) compute bpsi for boundaries
    // Note: including cornerpoints!
    // hence no concurrent kernels allowed with these kernels!
    common::setGpu1DVector<REAL>(m_dpsi_in_N_ext, hpsi_in_N_ext);
    common::setGpu1DVector<REAL>(m_dpsi_in_E_ext, hpsi_in_E_ext);
    common::setGpu1DVector<REAL>(m_dpsi_in_S_ext, hpsi_in_S_ext);
    common::setGpu1DVector<REAL>(m_dpsi_in_W_ext, hpsi_in_W_ext);

    // North
    if (m_bcN == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dpsi_in_N_cen, hpsi_in_N_cen);

        wavescomputer::solvepsi::reflectiveneumann::north<REAL>
            (m_dbpsi, m_dstoN, 
             m_dpsi_in_N_ext, m_dpsi_in_N_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcN == 1) // Sommerfeld boundary
    {
        common::setGpuBdVector<REAL>(m_dbufc, hpsiN);
        wavescomputer::reorder<REAL>(m_dpsiN, m_dbufc, m_Nx1);

        wavescomputer::solvepsi::sommerfeld::north<REAL>
            (m_dbpsi, m_dstoN,
             m_dpsiN, m_dpsi_in_N_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "solvepsi<north>: unknown boundary condition");
    }


    // East
    if (m_bcE == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dpsi_in_E_cen, hpsi_in_E_cen);

        wavescomputer::solvepsi::reflectiveneumann::east<REAL>
            (m_dbpsi, m_dstoE, 
             m_dpsi_in_E_ext, m_dpsi_in_E_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcE == 1) // Sommerfeld boundary
    {
        common::setGpuBdVector<REAL>(m_dbufr, hpsiE);

        wavescomputer::reorder<REAL>(m_dpsiE, m_dbufr, m_Nx2);

        wavescomputer::solvepsi::sommerfeld::east<REAL>
            (m_dbpsi, m_dstoE,
             m_dpsiE, m_dpsi_in_E_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "solvepsi<east>: unknown boundary condition");
    }


    // South
    if (m_bcS == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dpsi_in_S_cen, hpsi_in_S_cen);

        wavescomputer::solvepsi::reflectiveneumann::south<REAL>
            (m_dbpsi, m_dstoS, 
             m_dpsi_in_S_ext, m_dpsi_in_S_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcS == 1) // Sommerfeld boundary
    {
        common::setGpuBdVector<REAL>(m_dbufc, hpsiS);
        wavescomputer::reorder<REAL>(m_dpsiS, m_dbufc, m_Nx1);

        wavescomputer::solvepsi::sommerfeld::south<REAL>
            (m_dbpsi, m_dstoS,
             m_dpsiS, m_dpsi_in_S_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "solvepsi<south>: unknown boundary condition");
    }


    // West
    if (m_bcW == 0)      // reflective Neumann boundary
    {
        common::setGpu1DVector<REAL>(m_dpsi_in_W_cen, hpsi_in_W_cen);

        wavescomputer::solvepsi::reflectiveneumann::west<REAL>
            (m_dbpsi, m_dstoW, 
             m_dpsi_in_W_ext, m_dpsi_in_W_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcW == 1) // Sommerfeld boundary
    {
        common::setGpuBdVector<REAL>(m_dbufr, hpsiW);
        wavescomputer::reorder<REAL>(m_dpsiW, m_dbufr, m_Nx2);

        wavescomputer::solvepsi::sommerfeld::west<REAL>
            (m_dbpsi, m_dstoW,
             m_dpsiW, m_dpsi_in_W_ext,
             m_nx[0], m_Nx1, m_Nx2);;
    }
    else
    {
        ASSERT(false, "solvepsi<west>: unknown boundary condition");
    }
}


void CudaWavesComputer::adjustPsi(const REAL dt_new,
                                  const REAL dt_old)
{
    REAL beta  = dt_new / dt_old;

    wavescomputer::adjustpsi<REAL>
        (m_dpsi, m_dpsim1, m_dpsim2,
         beta, 
         m_nx[0], m_ny[0],
         m_Nx1, m_Nx2);
}


void CudaWavesComputer::shiftPsi()
{
    wavescomputer::shiftpsi<REAL>(m_dpsim2, m_dpsim1, m_dpsi,
                          m_nx[0], m_ny[0]);
}


void CudaWavesComputer::boundaryPsi(Array2D<REAL> *hpsi,
                                    const Array2D<REAL> *hpsiN,
                                    const Array2D<REAL> *hpsiE,
                                    const Array2D<REAL> *hpsiS,
                                    const Array2D<REAL> *hpsiW)
{
    // North
    if (m_bcN == 0)      // reflective Neumann boundary
    {
        wavescomputer::boundarypsi::reflectiveneumann::north<REAL>
            (m_dpsi, m_dpsiN, 
             m_dpsi_in_N_ext, m_dpsi_in_N_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcN == 1) // weakly-reflective boundary
    {
        wavescomputer::boundarypsi::weaklyreflective::north<REAL>
            (m_dpsi, m_dpsiN, 
             m_dpsi_in_N_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundarypsi<north>: unknown boundary condition");
    }


    // East
    if (m_bcE == 0)      // reflective Neumann boundary
    {
        wavescomputer::boundarypsi::reflectiveneumann::east<REAL>
            (m_dpsi, m_dpsiE, 
             m_dpsi_in_E_ext, m_dpsi_in_E_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcE == 1) // weakly-reflective boundary
    {
        wavescomputer::boundarypsi::weaklyreflective::east<REAL>
            (m_dpsi, m_dpsiE, 
             m_dpsi_in_E_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundarypsi<east>: unknown boundary condition");
    }


    // South
    if (m_bcS == 0)      // reflective Neumann boundary
    {
        wavescomputer::boundarypsi::reflectiveneumann::south<REAL>
            (m_dpsi, m_dpsiS, 
             m_dpsi_in_S_ext, m_dpsi_in_S_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcS == 1) // weakly-reflective boundary
    {
        wavescomputer::boundarypsi::weaklyreflective::south<REAL>
            (m_dpsi, m_dpsiS, 
             m_dpsi_in_S_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundarypsi<south>: unknown boundary condition");
    }


    // West
    if (m_bcW == 0)      // reflective Neumann boundary
    {
        wavescomputer::boundarypsi::reflectiveneumann::west<REAL>
            (m_dpsi, m_dpsiW, 
             m_dpsi_in_W_ext, m_dpsi_in_W_cen,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else if (m_bcW == 1) // weakly-reflective boundary
    {
        wavescomputer::boundarypsi::weaklyreflective::west<REAL>
            (m_dpsi, m_dpsiW, 
             m_dpsi_in_W_ext,
             m_nx[0], m_Nx1, m_Nx2);
    }
    else
    {
        ASSERT(false, "boundarypsi<west>: unknown boundary condition");
    }

    common::getGpu2DVector<REAL>(hpsi, m_dpsi, m_nx[0]); // needed for weakly-reflective
}


void CudaWavesComputer::timeDerivatives(const REAL hpslow,
                                        const Array2D<REAL> *hpressurePulse)
{
    common::setGpu2DVector<REAL>(m_dpressurePulse, hpressurePulse, m_nx[0]);

    wavescomputer::timederivatives<REAL>
        (m_ddt_zeta, m_ddt_phi,
         m_dU1, m_dU2, m_dh0, m_dD0,
         m_dzeta, m_dphi, m_dpsi,
         m_dSRCz, m_dpressurePulse,
         hpslow,
         m_Nx1, m_Nx2,
         m_dx1,m_dx2,
         m_cx[0], m_cy[0],
         m_nx[0],
         m_g);
}


void CudaWavesComputer::smoothZetaPhi(const REAL eps)
{
    wavescomputer::smoothzetaphi<REAL>
        (m_ddt_phi, m_ddt_zeta,
         m_ddt_phim1, m_ddt_zetam1,
         eps,
         m_nx[0], m_ny[0]);
}


void CudaWavesComputer::leapFrog(Array2D<REAL> *hzeta,
                                 const REAL dt_new,
                                 const REAL dt_old)
{
    wavescomputer::leapfrog<REAL>
        (m_dphi, m_dzeta,
         m_ddt_phi, m_ddt_zeta,
         m_dzetam1, m_dzetam2,
         m_dphim1, m_dphim2,
         m_dalfa_brk,
         dt_new, dt_old,
         m_nx[0], m_ny[0],
         m_Nx1, m_Nx2);

    common::getGpu2DVector<REAL>(hzeta, m_dzeta, m_nx[0]); // needed for maxwaves
}


/// returns device pointer to the RHS b
REAL* CudaWavesComputer::getDevicePtr_bpsi()
{
    return m_dbpsi;
}


/// returns device pointer to the solution psi
REAL* CudaWavesComputer::getDevicePtr_psi()
{
    return m_dpsi;
}



